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Abstract 

Numerous animals live in, and locomote through, subsea soils. To move in a medium dominated by fric- 
tional interactions, many of these animals have adopted unique burrowing strategies. This paper presents a 
burrowing model inspired by the Atlantic razor clam (Ensis directus), which uses deformations of its body 
to cyclically loosen and re -pack the surrounding soil in order to locally manipulate burrowing drag. The 
model reveals how an anisotropic body - composed of a cylinder and sphere varying sinusoidally in size 
and relative displacement - achieves unidirectional motion through a medium with variable frictional prop- 
erties. This net displacement is attained even though the body kinematics are reciprocal and inertia of both 
the model organism and the surrounding medium are negligible. Our results indicate that body aspect ratio 
has a strong effect on burrowing velocity and efficiency, with a well-defined maximum for given kinematics 
and soil material properties. 
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1. Introduction 

There are many examples of animals that live in 
particulate substrates which have evolved unique 
locomotion schemes [1J. Two common strategies 
observed in biological systems are an undulatory, 
snake-like motion El [3l SI and a "two-anchor" 
system H El El [El 03] CGI- An example of 
the former is the sandfish lizard which wiggles its 
body from side to side in order to effectively swim 
through sand 0. Similarly, smaller organisms like 
C. elegans have been observed to move quite effi- 
ciently via an undulatory motion through granular 
media (3j S) . In contrast, soft-bodied organisms 
that live in particulate substrates saturated with a 
pore liquid generally use a two-anchor system to 
burrow. In this strategy, one section of the ani- 
mal expands to form a terminal anchor, while an- 
other section of the animal contracts to reduce drag. 
Once the contracted section is conveyed forward in 



the burrow, it is expanded to form the next termi- 
nal anchor and the previous terminal anchor is con- 
tracted and shifted forward. 

The burrowing model presented in this paper is 
inspired by the two-anchor locomotion scheme and 
body geometry of the Atlantic razor clam (Ensis di- 
rectus). Ensis is comprised of a long, slender set 
of valves (i.e. the two halves of the shell) which 
are hinged on an axis oriented longitudinally to the 
animal, and a dexterous soft foot which resides at 
the base of the valves. The burrowing cycle of En- 
sis is depicted in Figure [T] (a). The animal starts 
with its foot fully extended below the valves (A). 
Next, it uses a series of four shell motions to make 
downward progress: (B) the foot extends to uplift 
the valves while the valve halves contract to force 
blood into the foot, inflating it to serve as a termi- 
nal anchor; (C) the foot muscles contract to pull 
the valves downwards; and (D) the valves expand 
in order to form a terminal anchor and begin the 
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cycle again. 

The uplift and contraction motion of the valves 
draw water towards the animal's body, unpack- 
ing and locally fluidizing the surrounding substrate 
ITT31 . The initiation of valve contraction causes lo- 
cal soil failure around the animal and the uplift ve- 
locity is on the order of the pore fluid velocity re- 
quired to induce a fluidized bed below the animal|j 
Although the animal is too weak to pull its shell 
through static soil (which exerts a resistance that 
linearly increases with depth lfT4lD to typical bur- 
row depths, fluidization dramatically reduces drag, 
resulting in resistance forces that are depth inde- 
pendent |[T3ll . The aim of this paper is to analyze 
the kinematic motion of the shell and demonstrate 
that reciprocal body deformations can produce uni- 
directional motion in a substrate of varying fric- 
tional properties. 

2. Model 

Figure [TJb) shows the geometry of the simpli- 
fied model organism and the dynamics inspired by 
Ensis. The body consists of two components: a 
long cylinder of length L and radius r(t), which ap- 
proximates the valves, and a sphere of radius R(t) 
attached to the cylinder, acting as the foot. The ra- 
dius of the cylinder, the radius of the sphere and 
the distance between the two are known functions 
of time dictated by the organism. The length of the 
shell, L, is considered constant. 

2.1. Kinematics 

Ensis' 's burrowing motions are often erratic; the 
animal may wait anywhere from a few seconds iTTTTl 
to many minutes between digging cycles. As a sim- 
plification, we consider periodic motion for both 
the foot and the shell. First, the radii of the cylin- 
der and the sphere are given by 

r - ao + a' cos(a>t) , R = a® - b' cos(cot) (1) 
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Figure 1 : Schematics of motion for (a) burrowing Ensis and 
(b) simplified model organism inspired by Ensis. 



where a) is the frequency of motion and ao is cho- 
sen as the mean radius of both the cylinder and the 
sphere. 

Another important factor in the kinematics of 
burrowing Ensis is the extension and retraction of 
the foot and its temporal relation to the movement 
of the valves. To model this motion, we impose 
a sinusoidally changing distance between cylinder 
and sphere that is out of phase with the expansions 
and contractions by n/2: 



d - do + d' sin(a»0- 



2.2. Volume conservation 



(2) 



'A full description of Ensis burrowing mechanics is be- 
yond the scope of this paper, but can be found in 1131 . 



In the live organism, the expansion of the foot 
is driven by fluid squeezed out of the shell. This 
fluid may be treated as incompressible at clam-like 
speeds and is contained in a closed loop. Hence, by 
conservation of volume, as the cylinder contracts, 
the sphere expands and vice versa. For small de- 



formations, the change in the total volume of the 
cylinder can be approximated as 



Vcyl 



InrrL InLcoaod sin(<x>f). 



(3) 



The sphere's volume varies with opposite phase 
and, again for small deformations is approximated 

as 



v S ph = Ana R ~ Anua^b' sin(wO- 



(4) 



Combining these two relations and applying con- 
servation of volume shows that the amplitude of 
the changing sphere radius is related to the defor- 
mations of the cylinder by b' = (L/2ao)a'. 

The initial mean void fraction of the surrounding 
medium, 60, is defined as the ratio of volume occu- 
pied by pore fluid to the total volume. The change 
in void fraction of the soil adjacent to the organism 
as the cylinder collapses relative to its mean value 
is given by 



(eo - e) cy i = 



nL(r 2 - <Zq) 

"Vcyl 

2nLa , 
— — a cos(a»?) 
' cyl 



(5) 



where "V is the characteristic volume of the per- 
turbed soil, the extent of which depends on the ge- 
ometry of the burrowing organism and initial soil 
properties [? ], and e is the instantaneous local 
void fraction. In the same manner, the void frac- 
tion change in the soil surrounding the sphere is 
given by 
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(6) 



2.3. Drag 



The drag force on an object moving through 
a saturated particulate medium depends (nonlin- 
early) on a number of parameters. In general, it 
can be written as 



F D = n(e)V a S 



(7) 



where p is a resistance coefficient which depends 
on local void fraction e, surface roughness, shape 
of the body, etc., V is the body's velocity, a is an 
exponent that varies with Reynolds number, and 5 
is a geometric parameter that is associated with the 
body's contact area with the substrate. 

For Newtonian fluids, expressions for fi in the 
limits of both high and low Reynolds number flows 
are well-known. The parameter a characterizes the 
velocity dependence in the drag expression. At 
high speeds (or high Reynolds numbers) the drag 
force is strongly dependent on velocity; however, 
as viscous effects increase and/or the substrate ex- 
hibits increasingly solid-like behavior, this veloc- 
ity dependence weakens. For inviscid Newtonian 
flows, a = 2, n is related to the dimensionless 
drag coefficient, Co, by /u = Cop/2 where p is 
the density of the fluid, and S is an effective cross- 
sectional area. At low Reynolds numbers, a - \ 
and ju is proportional to the dynamic viscosity of 
the fluid, fif where the constant of proportionality 
depends on the geometry of the body. For a sphere, 
li = 6nfif and S is the radius of the sphere. 

It should be noted that the analysis presented 
in this paper is predicated on defining fif as an 
effective viscosity which correlates shear stresses 
to strain rate. This behavior is markedly different 
than critical state granular shear flow, where inter- 
particle frictional interactions induce shear stresses 
that depend on confining pressure, and are rela- 
tively independent of strain rate lTT4l . Visualization 
of the fluidized substrate around burrowing Ensis 
show void fraction ranges of 0.42 < e < 0.46 |[T3l . 
which is unpacked beyond the point of incipient 
fluidization (e « 0.41) for uniform spheres |[22l . 
As such, the substrate surrounding burrowing En- 
sis can be modeled with a fluid-like viscosity that is 
a function of void fraction, for which there are nu- 
merous empirically-derived expressions in the lit- 
erature IT7HI81 [1911201. 

Following the sequence of events depicted in 
Figure [jjb), inward (collapsing) motion of the 
cylinder increases the local void fraction and con- 
sequently decreases the resistive drag force on the 
cylindrical body. In contrast, the drag on the sphere 



(which expands as the cylinder collapses) increases 
as the cylinder collapses. Assuming small local 
changes in void fraction Ae = e - £o about an initial 
eo, we can write the local resistance coefficient of 
the cylinder, jj. cy i{e) as 



Hcyiie) 
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where /j.' ; can be written as a function of the ge- 
ometry of the organism and the material proper- 
ties of the substrate by substituting equation ([5]) 
into d8|). In general, fi' l is assumed to be less 
than jUo otherwise the resultant negative drag on the 
moving body is not physical. For the sphere, fi sp h 
fluctuates out of phase with /i c> ,/ and is given by 
Msph = y"0 + H' sph cos(ajt). 

2.4. Burrowing velocity 

Since we are considering an inertialess limit, the 
total force on the model organism equals zero at 
every instant in time. This condition enables us 
to calculate a net burrowing velocity of the cylin- 
der/sphere system. 

The only forces acting on the system are drag on 
the cylinder and drag on the sphere which must be 
equal and opposite: 



HcylWcylV S cy l - fi sp h\V sp h\ a S sp h . 



(9) 



Since the distance between the cylinder and the 
sphere, d, is prescribed by the digger, the two ve- 
locities must be related via a kinematic constraint: 
d - V cy i - V sp h, where the dot indicates a time 
derivative. 

Solving for the velocity of the cylinder (which 
represents the shell of the digging clam) and substi- 
tuting the kinematics defined in ([2]) we find the di- 
mensionless digging velocity of the sphere, V sp h = 

v sph /(cod'y. 

COS(dtf) 



V S ph - ~ 
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Figure 2: Normalized burrowing velocity versus aspect ratio 
ao/L with different normalized soil parameters yu'/y"o- Dashed 
lines indicate a = 1 (strong dependence on velocity) and solid 
lines correspond to a = 0.1 (weak dependence on velocity). 



Figure [2] shows the relationship between this di- 
mensionless burrowing velocity averaged over one 
cycle, V, and the "shell" aspect ratio ao/L with 
varying normalized perturbed resistance coefficient 
n' /Ho and varying velocity dependence a. For sim- 
plicity, we have approximated {i' s h - //. = //, 

and S C yi and S sp h are chosen as 2naoL and 4na^, 
respectively. Note that as an approximation, choos- 
ing // h - fi' l is a reasonable first estimate but, in 
reality n' s h and fi' c tl are dictated by the constitu- 
tive equation relating // and e (which can be quite 
complicated and is often determined empirically 
for soils) and the geometry of the digger. 

Figure ^indicates that, for our chosen sinusoidal 
kinematics, the maximum burrowing velocity oc- 
curs at an aspect ratio of ao/L = 0.5 regardless 
of the material properties of the soil. As the nor- 
malized perturbed resistance coefficient fi' /hq is in- 
creased, the maximum velocity increases and, for 
small a, the velocity profile flattens out. These 
trends indicate that, as the resistance becomes more 
sensitive to changes in void fraction the burrowing 
velocity increases. Figure [3] shows a contour plot 
of burrowing velocity as a function of a normal- 
ized perturbed resistance coefficient // //jq and a at 
the fastest aspect ratio ao/L = 0.5. The burrowing 
velocity increases with decreasing dependence of 




Figure 3: A contour plot of maximum normalized velocity V cy i 
time averaged over one cycle as a function yu'/y"o and a at the 
optimal aspect ratio, a /L = 0.5. 



body forces on local velocities and increasing per- 
turbed resistance coefficient. 

2.5. Efficiency 

An actively burrowing animal consumes power 
as it deforms the surrounding medium and a certain 
fraction of that power is transformed into useful 
unidirectional motion. A typical hydrodynamical 
efficiency can be defined as the ratio of the power 
required to drag the digger through the soil at the 
average digging velocity (namely the useful frac- 
tion of the power) to the total power required to 
deform the substrate which can be expressed as 



fJ-oV (S C yl + S sp h) 
2ai=sph,cyl r £) ' 'i + "i 



(12) 



where overline indicates a time-averaged quantity 
and the latter terms in the denominator, P cy i and 
P sp h, indicate the power dissipated by an expand- 
ing (or shrinking) cylinder and sphere respectively. 
To evaluate the power associated with expansion, 
we first calculate the corresponding stresses as 
cr S ph = 2/UfdR/dr = (4/3)fifV sph /v sph and cr C yi = 
2nfdr/dr - fifV C yi/v cy i where v is the volume of 
the object and v is a dilation rate. As before, {if is 
the dynamic viscosity of the surrounding medium 
which is proportional to our resistance coefficient 



at low Reynolds numbers. The power dissipated is 
then given by crv or 
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(13) 



Using eqs. ([3]j4j> and approximating the normalized 
perturbed soil parameter as the same as the normal- 
ized perturbed radius {a' jao = ////io), we can eval- 
uate the power dissipated in expansion and contrac- 
tion and compute the efficiency of digging. Note 
that this approximation only considers the dissipa- 
tion associated with viscous stresses generated by 
an expanding object and at high Reynolds numbers 
there are additional components associated with Pj. 

Figure[4]shows the efficiency as a function of as- 
pect ratio ao/L. When either the cylinder or the 
sphere undergoes a large variation in size, the lo- 
cal viscosity fluctuates significantly, increasing the 
burrowing velocity. However, these large deforma- 
tions also dissipate a considerable amount of en- 
ergy and hence, in contrast to the velocity, the ef- 
ficiency decreases as the normalized perturbed soil 
parameter increases. In addition, larger normalized 
perturbed resistance coefficients broaden the effi- 
ciency profile as a function of ao/L. Figure[5]shows 
a contour plot of the efficiency as a function of a 
and // Ihq at the optimal aspect ratio, ao/L = 0.5. 
At small a, the efficiency is large owing to the 
corresponding high burrowing velocity. This high 
efficiency rapidly drops as a increases and varies 
weakly in the normalized perturbed resistance co- 
efficient. 

To compare this with other animals locomot- 
ing in a fluid environment, the hydrodynamic ef- 
ficiency is roughly 0.01 for small micro-organisms 
in a viscous fluid, and about 0.5 for large swim- 
ming animals in an unbounded fluid. Our re- 
sults indicate that burrowing animals have rela- 
tively high efficiencies despite the large resistivity 
of their surrounding environments. 
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Figure 4: Efficiency versus aspect ratio ao/L with varying nor- 
malized soil parameters yu'/yUo- Solid lines indicate a = 0.1 
(weak velocity dependence) and dashed lines correspond to 
a = 1 (strong velocity dependence). 



3. Discussion 

Given this formulation, there are a number of 
limiting cases that can be addressed analytically, 
yielding further insight into optimal geometries for 
burrowing. 

3.1. Limiting case 1: Small resistance perturba- 
tions 

In the limit /////o ^ 1> we can estimate y 1 ^* to 
first order in ///j"o as 
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were we have again approximated ju' c q ~ /j.' s h = 
//. Combining this with equation ( 10), we find the 
dimensionless instantaneous digging velocity can 
be represented as 
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where /(ao/L) = 1 + \S sp h/S C yi) ■ The mean di- 
mensionless digging velocity of the clam, namely 
V sp h time-averaged over one cycle, is given as 
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Figure 5: A contour plot of efficiency 77 as a function // //i 
and a with Oq/Z, = 0.5. 



Since / > 1, the body burrows downward in the 
vertical direction. 

In order to maximize digging velocity for a given 
geometry we set, 







dV 



f-2 



(16) 



d(S S ph/Scyl) P 

indicating that the maximum burrowing velocity 
occurs when / = 2 or equivalently when S sp h = 
S C yi regardless of the value of a or //. Thus, with 
S cyi = InaoL and S sp f, = 4nal the maximum ve- 
locity occurs at ao/L = 0.5. However, this optimal 
aspect ratio depends on our choice of S , which de- 
pends on the details of the burrowing system. 

3.2. Limiting case 2: Low Reynolds number New- 
tonian flows 

As the Reynolds number approaches zero, forces 
acting on the bodies are linearly proportional to ve- 
locity {a - 1). Expressions for these forces can be 
derived analytically, in particular, for an infinitely 
long cylinder aligned with the flow, the drag force 
is given by 2nfifVL and the drag force on a sphere 
is given by 67TjJ.fVao. We again consider the limit 
of small resistance perturbations // due to local 
changes in particle packing fraction in a viscous 
fluid. 

The calculation is the same as the previous cal- 
culation with the exception that S sp h = «o and 



S cy i - L/3. Hence the corresponding burrowing 
velocity is maximized at ao/L = 1/3, correspond- 
ing to a more elongated cylinder than in the previ- 
ous analysis^ 

3.3. Near random close packing 

For simplicity, our model has been linearized as- 
suming small variations in the packing fraction. 
This reduced model provides a number of gen- 
eral insights into how reciprocal motion of non- 
symmetric bodies can generate unidirectional mo- 
tion in a saturated soil. However, as the packing 
fraction of the substrate approaches critical transi- 
tional values (i.e. approaching random close pack- 
ing), a fully nonlinear viscosity model is more ap- 
propriate. 

There are a number of effective viscosity models 
which correspond to the special case of a - 1 as- 
sociated with burrowing at low Reynolds numbers. 
A few such models relating the effective viscosity 
to the packing fraction, <p, are itemized in the table 
below. 



V-lflf 


Reference 


1 +2.50 


Einstein [17] 


1 + 2.50 + 7.6<£ 2 


Batchelor & Green HT8TI 


(9\ WM L/i 
U) i-(<t>/<t> m ) i/3 


Frankel & Acrivos lfl9l 


\ <t>m) 


Krieger & Dougherty (20l 



2 Note that there is a subtlety in this calculation that needs 
to be addressed. The previous calculation, which should also 
be relevant at low Reynolds numbers, yielded an optimal as- 
pect ratio of 1/2, not 1/3. To rationalize this apparent discrep- 
ancy, consider the drag on a sphere in a low Reynolds number 
flow. There are two common (and equivalent) ways to express 
the drag force: F D = 6npfVa or F D = C D p/2V 2 S where 
Co = 24 1 Re and S = na^. In the first case, in our formulation, 
we set p = 6nfif and S sp i, = a Q resulting in an optimal aspect 
ratio of 1/3. In the second, we set/j = C D p/2 and S sph = S 
yielding an optimal aspect ratio of 1/2. To determine which is 
correct, we need to consider the constitutive relationship for 
the resistance coefficient. If 6npf can be well-approximated 
as a linear function of e, the first formulation is relevant. If, 
on the other hand C D p/2 = 6pf/(VR) (which includes the ge- 
ometric parameter R which also affects the void fraction) can 
be considered a linear function of e, the second estimate is 
more appropriate. 



If we repeat the previous calculation to determine 
average digging velocities, this time taking into ac- 
count the full nonlinear dependence cited in EOl 
and using typical parameters (rj = 2.5, <p m = 0.67), 
the resultant velocity becomes V = 0.1295, which 
is commensurate with the values predicted by the 
linearized theory. 

3.4. Previous numerical results 

While to the best of our knowledge this paper 
represents the first theoretical analysis of a simple 
burrower using local fluidization to propel itself, 
this type of digging strategy has previously been 
studied numerically by Shimada et. al. lfl6l . In 
that study, the authors used an event-driven gran- 
ular simulation to model a "pushme-pullyou" con- 
sisting of two expanding and contracting disks sep- 
arated by a spring. Both halves of the body were 
disks, hence aspect ratio was not a parameter in 
their study. In both studies (present and previous 
numerics) the burrowing velocity was found to be 
proportional to u at low frequencies. In the event- 
driven simulations, the authors found that this re- 
lation peaks at a critical frequency and the veloc- 
ity declines beyond this critical value. Our cur- 
rent theory is unable to predict these nonlineari- 
ties observed at high a) on account of the assump- 
tions made in the constitutive relationships, namely 
ju oc e. A more realistic constitutive model (which 
would depend on the details of soil type, prepara- 
tion, etc.) is likely to exhibit behavior that is quali- 
tatively similar to the numerical simulations at high 
frequencies. 

3.5. Conclusions 

In this paper, we introduce a simple theoretical 
model to capture key physical aspects of burrowing 
Ensis and other biological or engineered burrow- 
ing systems. Even though the cylinder and sphere 
motion is actuated reciprocally in an over-damped 
environment, net unidirectional motion is achieved 
because of varying drag on the bodies owing to 
local changes in void fraction. We find that bur- 
rowing velocities depend on the aspect ratio, ao/L, 
and that the "best" aspect ratio (i.e. the one that 
maximizes the velocity or the efficiency) depends 



on the geometric details of the drag force expres- 
sion. It is interesting to note that, while we found 
an optimal ratio of 1/3 for viscously dominated 
substrates, live razor clams have an aspect ratio 
closer to oq/L ~ 1/6. This discrepancy is likely 
to arise due to an over-simplification of the consti- 
tutive relationships describing the substrate or may 
be an indication that razor clams have not evolved 
to maximize digging speeds. To better describe the 
dynamics of Ensis, future work will focus on more 
precisely determining the parameters and e depen- 
dence in the force relation expressed in equation 
Q and investigating other cost functions. 
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